function [xh, Sx, pNoise] = ukf_measurement(state, stateCov, pNoise, U1, modelObject)

Xdim  = modelObject.statedim;                                % extract state dimension
Vdim  = modelObject.Vdim;                                    % extract process noise dimension

% Get and calculate CDKF scaling parameters and sigma point weights
h = modelObject.scaleFactor;
hh = h^2;

W1 = [(hh - Xdim - Vdim)/hh   1/(2*hh);                  % sigma-point weights set 1
    1/(2*h)                 sqrt(hh-1)/(2*hh)];

Zeros_Xdim_X_Vdim = zeros(Xdim,Vdim);
Zeros_Vdim_X_Xdim = zeros(Vdim,Xdim);

nsp1   = 2*(Xdim+Vdim) + 1;          % number of sigma points (first set)

Sx = stateCov;             % matrix square root of state covariance
Sv = pNoise.cov;         % matrix square root of process noise covariance
%------------------------------------------------------
% Sigma points generation

%Z is the sigma points,
Z   = cvecrep([state; pNoise.mu],nsp1);

%Sz is the augmented coveriances
%Sz =  [Sx 0
%      0  Sv];
Sz  = [Sx Zeros_Xdim_X_Vdim; Zeros_Vdim_X_Xdim Sv];
%hSz is the scaled sigma point set

%scale Sz by the scale factor
hSz = h*Sz;
%build the sigma point set (3.225 in the paper), Z is the sigma point
%set the sigma points
hSzM = [hSz -hSz];

%First state + the augmented sigma points
Z(:,2:nsp1) = Z(:,2:nsp1) + hSzM;


%------------------------------------------------------
% TIME UPDATE

%propagate sigma points through our time update equation (3.226), UU1
%is our input vector
X_ = f_fun( modelObject, Z(1:Xdim,:), Z(Xdim+1:Xdim+Vdim,:), U1);  % propagate sigma-points through process model

%Average our propagated sigma points with our weights (3.227)
xh_ = W1(1,1)*X_(:,1) + W1(1,2)*sum(X_(:,2:nsp1),2);

%set the state
xh = xh_;

%Sx, new covariance matrix, (3.228)
% start+1 :mid - mid+1-end
A = W1(2,1) * ( X_(:,2:Xdim+Vdim+1) - X_(:,Xdim+Vdim+2:nsp1) ) ;
% start+1 :mid - mid+1-end - 2 * first element
B = W1(2,2) * ( X_(:,2:Xdim+Vdim+1) + X_(:,Xdim+Vdim+2:nsp1) - cvecrep(2*X_(:,1),Xdim+Vdim));

%augment both coveriances and take QR dcomp.
[~,Sx_] = qr([A B]',0);

%set the covariance
Sx= Sx_';



